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ABSTRACT 



A full-wave solution has been obtained for radio wave propagation in the 
presence of an elevated tropospheric layer. This analysis was performed 
as part of the CNI (Communication, Navigation, Identification) experimental 
program. The tropospheric layer is modeled as a trilinear refractive 
index profile with a sufficient lapse rate so as to result in an elevated duct. 
The analytical solutions give the received signal level for air-to-air propa- 
gation paths at UHF, the corresponding signal fading level, and the space 
diversity distance to insure good quality reception. The analysis performed 
here can be applied to ionospheric propagation and to the underwater 
acoustic channel. 
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SECTION 1 



INTRODUCTION 

1- 1 Motivation for the Problem 

In the course of designing a propagation experiment for the 

■k 

Communication, Navigation and Identification (CNI) Program [1] , for 
air-to-air and air-to-ground communications links at UHF, it was recog- 
nized that anomalous tropospheric conditions could lead to serious de- 
gradation (or enhancement) in signal strength and to multipath (inter- 
ference) effects similar to those produced by the earth. These facts 

prompted an initial investigation into the subject of tropospheric 

** 

propagation. It was learned that tropospheric layers in the 

first few kilometers above ground are relatively common, but their 
precise effect on signal strength, fading characteristics, Doppler 
spectrum, available bandwidth, space diversity and relative terminal 
heights has been investigated in relatively few cases. 

Propagation in an anomalous atmosphere has been studied exten- 
sively using the theory of geometrical optics. The theory was quite 
successful in explaining the gross features of the results of air-to- 
air propagation tests at UHF [2]. However, the theory is approximate 
and it suffers from the drawback of all geometrical optics treatments 
in that the solution breaks down at the caustics, and furthermore, 
cannot predict field strengths beyond the horizon. Analytical solu- 
tions or full wave solutions for elevated layer propagation have re- 
ceived only a few approximate and incomplete treatments. Furry [5] has 
obtained approximate solutions which are valid for the trapped or ducted 
waves only, Northover [6] has treated (at VHF) the elevated layer as a dis- 

* A bibliography is listed at the end of the report. 

** The term layer means a stratum, of vertical thickness from a few 
meters to several hundred meters, within which the mean vertical 
gradient and/or the variance of refractive index are much greater 
than elsewhere (see Section 1.2). 
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continuity of the index of refraction, which is an unrealistic model 
in the general case, and Chang [5] considered (at VHF) only the trapped 
waves, essentially ignoring the effects of the earth. 

As a result of this initial investigation, it was seen that the 
present insufficient and incomplete knowledge of wave propagation at 
UHF in the presence of elevated layers could not lead to the effective 
planning of a propagation experiment. Moreover, only more precise 
analytical results along with a more realistic model of the layer would 
lead to a better interpretation of the final test data from such an 
experiment. It is with this objective in mind, that the present work 
was undertaken. 

1. 2 A Discussion of the Problem 

We will give in this section an outline of the problem and the 
results we hope to obtain. We will discuss the assumptions that we 
are making, and delineate the part of the work that is original. We 
will also include a description of the troposphere and the effects of 
elevated layers on received signal strength. 

The problem considered is that of determining the electromagnetic 
field produced by a vertical magnetic dipole (transmitter) in a radially 
inhomogeneous medium (troposphere) on or above a spherical earth. A 
full wave solution is obtained for the case when the variation of the 
index of refraction of the atmosphere is such as to produce an elevated 
tropospheric duct. The solutions for the field intensity are obtained 
as a function of frequency, layer height and vertical thickness, lapse 
rate or change of index of refraction with height (dN/dh) within the 
layer, and as function of transmitter and receiver heights and separa- 
tions. In the solution of the problem the WKB or phase integral solu- 
tions for a second order ordinary differential equation are used. The 
vork of previous authors is extended in three important ways. First, we 
use the "extended" WKB solutions of Langer [6) . This approach allows 
tne closing of the "gap" in the regular WKB solutions. It involves Hankel 



functions of order one-third. Second, the Langer solutions valid only 
for a single transition point are extended to be valid for three tran- 
sition points. This is essential to obtain results valid in the elevated 
duct and it is accomplished by obtaining solutions in three regions and 
joining these at the bottom and top of the elevated layer heights. And 
third, the resultant solutions for the trapped and leaky waves, are ob- 
tained from the solution of a single transcendental equation. This re- 
sult is, to the best of the author's knowledge, entirely original. 

The refractive index profile of the troposphere is given by three 

linear line segments of different slope. While at first sight this choice 

might seem artificial, it has been shown by ref ractometer measurements that 

index of refraction profiles do take on such abrupt changes and, hence, 

the trilinear profile becomes a valid approximation. In reality, all that 

is required for the presented solutions to apply is that the profile be 

"sufficiently" linear in the three segments . The second assumption 

we are making is to replace the spherical earth by a flat earth and modify 

the index of refraction, so called earth flattening approximation. In 

this case the solutions are valid to within two percent up to ranges of 

about half the radius of the earth, and this independent of the frequency. 

On the other hand, the error in the height-gain functions increases with 
5/2 

elevation as h and is proportional to the frequency [7] . At a 
frequency of 300 MHz and an elevation of 5000 feet, for example, the 
height-gain function is only in error by less than two percent. The 
third assumption is that the earth is a perfect conductor. This 
is a good approximation for horizontally polarized waves and for the low 
elevation angles of interest here. 

We will now turn to a discussion of the troposphere and the effects 
of elevated ducts on communication links. 



The sufficiency condition is given in Section 1.3. 
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The elevated inversion layer or elevated duct is characterized by 
a rapid decrease in the index of refraction, n over a relatively short 
height interval. The index of refraction n can be derived from basic 
thermodynamic laws and can be shown to give 



n = 1 + 



77.6 , 4810e. 
— (P + — m — ) 



10 



-6 



where p is the atmospheric pressure (millibars), T is the temperature 
(absolute scale) and e is the partial pressure of water vapor (milli- 
bars). The refractive index of the troposphere does not to a first 
order approximation depend on the frequency for wavelengths longer than 
1 cm. For waves of the millimeter range considerable losses appear, 
which can be considered by means of the introduction of a complex 
dielectric constant of air. Then the index of refraction n will depend 
on frequency. In this report we are not interested in frequencies higher 
than 10 GHz. 

In view of the small departure of n from unity (on the order of 

-4 6 
10 ) it is often convenient to use the quantity N = (n - 1) 10 , called 

ref ractivity . Near the earth's surface the index depends on the climatic 

and meteorological conditions and varies within 260-460 N-units. On 

the average, the value of N changes linearly with height, and for the 

middle latitudes the gradient of the change N with height h is 

^ ~ -AO N-units /km . 
an 

An elevated layer is shown qualitatively in Figure 1.1 at an altitude 
between 1000 and 1100 meters. 
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1000 




N-units 



360 



FIGURE 1.1 ELEVATED DUCT — 



N-PROFILE 



Most frequently these layers are caused by a mass of warm dry air over- 
lying a cool moist air mass. The usual occurrence of this phenomenon 
is at heights of approximately 1 to 2 kilometers. When the refractivity 
gradient exceeds 157 N-units/km, within the layer, the layer 

becomes an elevated duct. The changes in N have been measured to be as 
3 

high as 5 x 10 N-units per km. Sometimes 20-25 N-unit changes have 
been measured across a layer having a thickness of a few meters only. 

Even though, the modified index or earth flattening procedure has 
certain limitations as pointed out above, it offers great simplification 
and a large insight into duct propagation. For this reason, we will use 
the concept of the modified index of refraction to give a qualitative 
description of propagation in elevated ducts. The modified index of 
refraction seeks to accomplish two aims: (1) its vertical gradient 
should vanish at any elevation for which the path of a ray launched 
horizontally is a circular arc concentric with the surface of the earth; 



and (2) it should be readily calculable from the refractivity N. 
For this reason the modified index is defined by [8], 



M(h) = N + 10 6 - 
a 



where h is the height above the earth and a is the radius of the earth, 



height 
(m) 



1100 + 
1000 




360 



M-units 



FIGURE 1.2 



ELEVATED DUCT 
M - PROFILE 
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In a distribution in which M is constant rays are straight, when dM/dh 
is negative rays curve downward and when dM/dh is positive, the rays 
are refracted upward. The net result is that through the use of M(h) 
a spherical coordinate problem is more simply analyzed in rectangular 
coordinates, where plane wave modes have replaced spherical-harmonic 
modes. An M-profile roughly corresponding to the N-profile of Figure 
1.1 is illustrated. We will use a profile similar to Figure 1.2 in the 
solution of our problem. 

The disturbing effect of ■ tropospheric ducts on air-to-air and air- 
to-ground communication links has been investigated in several important 
experiments [9, 10, 111. In general three types of disturbances 

have been observed which cannot be explained in terms of ground 
reflections and which are evidently associated with the tropospheric 
propagation medium. These disturbances are (i) intervals of range in 
which average signal has decreased 15 or 20 dB (radio-holes) , (ii) 
intervals of range in which the average signal exceeds the free-space 
value by 10 or 12 dB (anti-holes) and (iii) intervals in range in 
which interference type of lobes (multipath) of large amplitude (10- 
20 dB) appear. It has been observed that radio holes and anti-holes 
occurring far from within to far beyond the standard horizon are very 
common air-to-air. For instance, all of the 50 air-to-air propaga- 
tion flights made by Wright Air Development Center in 1951-1953 en- 
countered radio holes; these flights were made under fair weather 
during May-October in a three year period in Ohio [9l . Radio holes 
far within the standard horizon are not so common ground-to-air. 
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1.3 Outline of Approach 

The procedure is to reduce the problem of determining the electric 
field strength of a dipole anywhere in space to the solution of an in- 
homogeneous Hertzian potential wave equation with radial and angular 
coordinate dependence. The solution is subject to the boundary condition 
at the earth, the Sommerfeld radiation condition, the continuity of the 
solution plus its derivative at all non-singular points and the condition 
that the solution must remain finite at the dipole source. By separation 
of variables the Hertzian potential is represented as an infinite sum of 
radial and angular dependent functions. The angular functions are 
Legendre functions and the radial functions are as yet unspecified. The 
coefficients of the summation can be obtained from the boundary conditions. 
The resulting series is, however, a very slowly converging series. To 
improve the rate of convergence of this series the Watson transformation 
is applied. The individual terms of this resulting series are evaluated 
at the roots for which the Wronskian vanishes. This will give a general 
solution. 

The refractive index profile is then chosen as to give three trans- 
ition points; this corresponds to an elevated duct. The radial dependent 
functions mentioned above are evaluated by an approximate solution to the 
radial wave equation with the "extended WKB" method due to Langer [6], 
at each transition point and in three separate regions. The solutions 
are Bessel functions. The solutions in the three regions are then joined 
by the continuity conditions. The fact is then emphasized that the 
propagation problem reduces to the solution of an eigenvalue problem for 
an ordinary differential equation. Hence, the solution has to satisfy 
the boundary condition at the surface of the earth and the radiation 
condition at infinity. The resulting set of complex eigenvalues and the 
associated eigen-f unctions represent a complete set of solutions to the 
problem. Substituting this set of eigenvalues and eigen-f unctions in 
the general series solution, we can determine the field strength due to 
the radiating dipole at any point in space. 
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In the solution of this problem it will be seen that an eigen- 
value ordinary differential of the form occurs: 

d u ~ 

+ k Z [y(h) + A ] u =0 (1.1) 
,,2 o n n v ' 

dn 

where A^ and u^ are the eigenvalues and eigenf unctions , respectively , 

k = 2tt/A, and y(h) is related to the modified index of refraction by 

y(h) = 2 x 10 (M(h) ~ M ) . The functions u must satisfy the boun- 

o n 

dary condition at the earth and the radiation condition at infinity. 

The form of the solution to this differential equation will depend upon 

the term in brackets and in particular upon the eigenvalues A^ which 

are in general complex. When the A^ have imaginary parts (positive) 

which are comparable in magnitude to their real parts, the corresponding 

wave is said to be "leaky". This is because the value of the function 

u^(h) is found to increase rapidly with height. When the imaginary 

part of A becomes small and the real part of A is comparable to 
n n 

(-M(h c )) (see Figure 1.1), we have "transitional" waves. And when 

the imaginary part of A^ ^ -M(h^), the waves are "trapped", because the 

functions u (h) has its largest value between h and h . We will use 
n a c 

the "extended" WKB (Wentzel, Kramers, Brillouin) solutions of equation 
(1.1) in the three cases of leaky, transitional, and trapped waves. 

In general, the term [y(h) + A ] vanishes at three points which, 
though complex, are extremely close to the real values h^, h^ , h^ shown 
in Figure 1.1. 

Such points are called turning points or transition points. This 
description has its origin in quantum mechanics where an equation very 
similar to equ. (1.1), namely Schrbdinger 1 s equation, has been the sub- 
ject of considerable analysis. The solutions of this equation make a 
transition from, say oscillatory to exponential at the transition point, 

The Jeffreys, Wentzel, Kramers, Brillouin (JWKB) or WKB method is 
a procedure for finding approximate solutions to a second order differ- 
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M 



FIGURE 1.3 



M-PROFILE 



ential equation of the form 



u" + f (z)u = 



(1.2) 



in some region, R, in which f(z) is analytic and contains no zeros, 
and in which f(z) varies sufficiently slowly. This means that the 
first two derivatives of f have to be sufficiently small throughout 
the region in question, with f ^ therein [4]. The WKB solution of 
(1.2) is given by 
1 



u = f 



(Ae 



iF(z) 



+ Be 



-iF(z) 



(1.3) 



where 




z 



(1.4) 



o 
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and A and B are arbitrary constants. However, the WKB approximation 

always fails sufficiently near a zero of f(z). A form of the solution 

of (1.2) which is valid up to the zero of f (z) itself was given by 

R. E. Langer [6]. This solution has the form 

1 

(1.5) 



u = (fr> 



A fl^ (F) + B (F) 



where F = F(z) is as defined above and A and B are arbitrary constants, 
and the are Hankel functions of the first and second kind. This 

approximation satisfies the equation 



u - + / Q 2 _ 1_ ^ + 3 g± 

V 36 F 2 4 Q 2 




- f n u = d- 6 ) 



2 

where Q = f . Hence, if we can show that 



1 2li _ 5_ £ _ 1 g» 2 
h Q 2 36 p 2 2 Q ^ 



(1.7) 



the solution in (1.5) is a very good approximation of the original 
differential equation (1.2). 

A limitation of the Langer WKB solution extension is that the 
method applies only to a single transition point. Hence, in our 
approach to the problem we shall obtain "extended" WKB solutions of 
the differential equation in a region including only one transition 
point. These solutions are then "joined" at the border of each region 
by employing the continuity of the functions and their derivatives » 
so that the general solution is valid in the semi-infinite interval 
extending from the surface of the earth to infinity. 



Northover [5 ], p. 37ff, shows that for the WKB solution to hold 
we must have If'/ f 3/2 T « i a nd f"/f 1/2 « 1. 
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1.4 Contents of the Report 

We will give here a brief outline of the remainder of this report. 
In Section 2, we formulate the problem in a precise mathematical way. 
The three dimensional vector solutions of Maxwell's equations are re- 
duced to a one dimensional vector by the introduction of the Hertz 
potential. By placing the transmitting source on the polar axis, ad- 
vantage is taken of the symmetry of the problem. As a result, the 
Hertz potential has only an r and 6 dependence. The electric and 
magnetic fields can then be written in terms of this Hertz potential. 
In Section 3, we obtain the slowly converging series solutions for the 
Hertz potential. By separation of variables this series can be written 
as the product of the radial dependent , r and angular dependent , 
functions. The angular functions are Legendre polynomials, and the 
radial functions are as yet unspecified; since they depend on the 
particular index of refraction profile assumed. In Section 4, we apply 
the Watson transformation to obtain a more rapidly converging solution. 
This is done by expressing the series solution for the potential as a 
contour integral in the complex plane enclosing an infinite number of 
singularities. By Cauchy's residue theorem, the more rapidly converging 
series for the potential is now a sum over the residues enclosed by 
the contour. In Section 5, we use Langer's "extended" WKB solutions 
to obtain approximate solutions of the radial dependent functions. And 
in Section 6, we obtain the radial dependent solutions for three trans- 
ition points. In Section 7, we transform these radially dependent sol- 
utions from that of a spherical geometry into a "flat earth" geometry. 
In Section 8, we present a model for the layered troposphere. And in 
Section 9, we obtain the expression for the Hertz potential and the 
electric field. In Section 10 we give the numerical results for the 
potential. The discussion in this section includes details of the numerical 
evaluations. Section 11 contains a summary. Relegated to the Appendices 
are the mathematical details of the contour integral representation and 
the modified Hankel functions. Also included there is the justification 

12 



for representing the three dimensional vector solutions of Maxwell's 
equations by the one dimensional Hertz potential. 
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SECTION 2 

PRECISE MATHEMATICAL FORMULATION OF THE PROBLEM 



We assume the earth to be a sphere of radius a, and introduce spher- 
ical polar coordinates r, 0,? as defined in Figure 2.1. We are given 
a vertical magnetic dipole of moment ml with a harmonic time var- 
iation e at a point D(b, 0, 0) and we are required to determine 

the electromagnetic field at an arbitrary point P(r, 0, 9) on or above 

2 2 1/2 
the surface of the earth and at a distance R = (r + b - 2br cos 0) 

from D. We assume the earth to be immersed in an inhomogeneous medium 

whose dielectric constant k is a function of radial distance r only, 

** 

k = k(r); all effects of the ionosphere have been neglected . Inside 

the earth k(r) is assumed to be a constant, k£. In any medium the 

propagation constant is i< = [iuy (o + iaje)]!/^ so that in air k = 2tt/X, 

1/2 ° 
while in the earth k 9 - k (K - 60 iaX) ' , where K is the relative 

L o 

dielectic constant, a is the conductivity (in mhos per meter) and the 

permeability y is assumed to be the same inside and outside the earth. 
2 2 

In this problem k (r) = oj ye(r) for the troposphere. 

The electromagnetic field must satisfy Maxwell's equations, both in- 
side and outside the earth, it must possess at the point D a singularity 
corresponding to that of the dipole, the electric field and magnetic 
field must be continuous everywhere, while at infinity the field vectors 
must satisfy the Sommerfeld radiation condition. 

Let E and H represent the electric and magnetic field vectors, 
respectively; the factor exp (io>t) being understood throughout. Then 
Maxwell's equations in m.k.s units, become the following: 



The material in this and the following sections is part of the author's 
dissertation for the Ph.D. at the University of Pennsylvania, 
Philadelphia. 

** It has been shown experimentally that at frequencies higher than 
30 MHz the ionosphere has negligible effects on terrestrial wave 
propagation. 
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V x E = -iojpH V • H = 

(2.1) 

V x H = iweE V • eE = 0. 

It has been shown by Friedman [ 12 ] that in spherical coordinates the 
electromagnetic field given by (2.1) can be obtained from a Hertz 
vector n which has only its radial component non-zero, even when the 
dielectric constant is a function of r. Since the dipole is located 
on the polar axis, symmetry exists, and the vector II will not depend 
on 9 but only on the coordinates r,6. We define II = (rU(r,9) ,0,0) . 
In terms of U the electric and magnetic fields are: 



(2.2) 



311 

V - E e =0 - V° \ " 36 

3r 

Equations (2.2) represent the field of a magnetic dipole and are sat- 
isfied for 

-ikR 

U-M^L_ (2.3) 

where M is proportional to the magnetic moment of the dipole, its 
value being 

M-~~ (2.4) 



where I is the electric current in the loop of area S and b is the 
distance from the center of the earth to the dipole. 
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The scalar wave function U(r,0) must satisfy several conditions. 
First U(r,0) must satisfy the partial differential equation 

I. V 2 U + k 2 U =0 (2.5) 

in the exterior region r > a, except at the dipole, and in the interior 
region r < a. Second, 

II. U and |^ (rU) (2.6) 

must be continuous everywhere except at the source. Thirdly, as R 

-lk (b)R 

This means that in the neighborhood of the dipole the potential U be- 
comes that of an oscillating dipole in an homogeneous medium whose 
permittivity is equal to the local value at the dipole position. 
Finally, U must satisfy the Sommerfeld radiation condition given by 

IV. lim r (— + ik U) = . (2.8) 

dv O 

The fact that U(r,6) possesses a singularity at the dipole, r = b, 
6=0, while otherwise it satisfies the wave equation (2.5), can be 
expressed very concisely using the Dirac 6-function. We write 



V 2 U+k 2 U=-M 6(r - b)6(9) (2.9) 
2Trr sin 

2 

where M is the source strength and the factor 27Tr sin is designed 
to normalize the solution. 
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SECTION 3 



SERIES SOLUTION OF THE SPHERICAL HERTZIAN WAVE EQUATION 

We shall now obtain a general solution to the inhomogeneous 
Hertzian potential wave equation given by (2.9), subject to the re- 
quired boundary conditions discussed in the previous section. Rather 
than writing the total field as the sum of a primary field due to 
the dipole source and a scattered field due to the earth, we shall use 
a self-consistent method. In particular, we will show that the solu- 
tion to the inhomogeneous wave equation (2.9) has the form of an ex- 
pansion in zonal harmonics, 

oo 

U(r,0) =]C i^h u n (r) P n (cos 0) (3.1) 
n=0 

where n are integer values and u n ( r ) are solutions to the equation 



d^(ru ) 

. 2 ( r ) _ n(n+l) 



H k2 



or 



2 - * \w 2 
dr L r 



To show this we use the orthogonality property of the Legendre poly- 
nomials P (cos 0) in (3.1) and find that 
n 



(ru n ) = ^~ 6 (r-b), (3.2) 



u n (r) = / U(r,6) P n 



(cos 6) sin 6 d6 . (3. A) 





If we multiply equation (2.9) (for spherical coordinate system) by 

P n (cos 6) sin 6 d6 and integrate between and ir, this equation transforms 

into the ordinary differential equation (3.2) • 
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The foregoing implies that for U(r,0) given by (3.1) to be a 
solution of the inhomogeneous wave equation (2.9), the coefficients 
u^(r) must satisfy eq. (3.2). Furthermore, since the Legendre poly- 
nomials form a complete set, U(r,0) as given by (3.1) is uniquely 
specified. Hence, U(r,0) given by eq. (3.1) is the total solution 
of the inhomogeneous wave equation. 

The boundary conditions that had been imposed on U in Section 2 
must now be satisfied by ( ru n ) > namely 

1. (ru ) and (ru ) must be continuous 

n ar n 

2. must be regular at r = 

3. The Sommerfeld radiation condition 



lim 1^- (ru ) - ik (ru )| = 
1 dr n on 1 



must be satisfied. 

In order to solve equation (3.2) we must know the solutions of 
the homogeneous differential equation 



^2 (rU n } + 

dr 2 



k 2 (r) _n(n±il 



r 



(ru ) = . (3.5) 
n 



2 2 

When k (r) = k Q , a constant (case of airless atmosphere), the solutions 

of (3.5) are linear combinations of the Hankel functions of half-integer 

(1) (2) 2 

order, namely H ,,/«(*) and H , /0 (r), If k (r) is not a constant, the 
n+l/z n+l/z 

solutions of (3.5) behave like Bessel functions, but their exact behavior 
will depend on the variation of k(r) with r. 

3.1 General Solution for the Radial Functions 



We let 



u (r) = a V (r) for r < b 
n n n — 

u (r) = b Q (r) for r > b 
n n n — 
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where the v n (r) and Q^*") are independent solutions of the homogeneous 
differential equation in (3.5). They, furthermore, satisfy the re- 
quired boundary condition in their respective domains. The general 
solution (3.1) is then given by 



00 

U(r,0) = /I C ^T i) p < cos 9 > 



n=0 



a V (r) 
n n 



b Q (r) 
n x n 



r < b 



(3.6) 



r > b . 



The coefficients a^ and b^ may be determined by applying the required 
boundary condition at r = b. At r = b, the functions (* u n ) are con- 
tinuous and their derivatives have a junp of magnitude -M/2tt b. From 
the definition of v n ( r ) an ^ above, the functions u n (0 are re- 

gular at r = and satisfies the radiation condition. Then we can 
prove that 



a = 



-M 

2ttW x n 
n 



<L0>> 



(3.7) 



b = 



-M 

2ttW "n 
n 



V (b) 



(3.8) 



where W is the Wronskian of (rQ ) and (rV ) evaluated at r = b. Sub- 
n x n n 

stituting (3.7) and (3.8) in (3.6), we obtain 

Q (b) V (r) r < b 
n n — 

(3 

.Q n (r) V n (b) r > b . 

This is the solution to our problem in its most general form for an 
arbitrary variation of the index of refraction of the atmosphere. 
Equation (3.9) represents Greens 1 function for the inhomogeneous wave 
equation (2.9). The general requirement of the reciprocity of Greens' 
function is satisfied owing to the fact that eq. (3.9) is symmetric in 
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oo 

n=0 n 



r and b. The reciprocity with respect to angles 8 and 0^ (we con- 
sidered the case in which the latter is zero) can be expressed by re- 
placing cos 8 by cos (8 - ) which is symmetric in and q . It has 
been shown that the solution (3.9) can be used for numerical computa- 
tion only if ka is small compared to unity. For example, at 100 MHz 
about a million terms are required to evaluate the series. In the 
following sections the solution (3.9) will be transformed so that it 
can be readily computed for large values of ka. 
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SECTION 4 



THE WATSON TRANSFORMATION 



The series expression for U(r,0) is very slowly converging. For 
example, for ka of the order of 10^, where k = 2ir/X and a is the radius of the 
earth, approximately 10^ terms of the series are required to evaluate 
the value of U(r,0) at one point. Such numerical calculations are 
even prohibitive in this age of high speed computers. A way around 
this laborious task was found by Watson [13]. Watson showed that the 
series (3.9) can be expressed as a contour integral in the complex plane 
enclosing an infinite number of singularities. By Cauchy's residue 
theorem this integral is then equal to the sum of the residues enclosed 
by the contour. The net effect of the Watson transformation is to have 
a much faster converging series. The difficulty that is left, is to 
locate the singularities, namely zeros of the Wronskian (see Appendix B) . 
From eq. (B.12) we have the series representation 



u(r.e) -f 22 



sin vir dW 
as 

s=v 



Q v (b) V v (r) r < b 

Q v(r) V^Cb) r >_ b 
(4.1) 



where v represents the zeros of the Wronskian W g and the derivative 
(d/ds) is with respect to the complex variable s. 

The residue 

v + \ P (cos Or - 0)) 

Q (b) V (r) (4.2) 



where 



sin vir _ _ v v 

w 

v 



W ' = ~- W | (4.3) 
v ds s 1 

s=v 
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can be interpreted in another way. Since the Wronskian W g = 0, the 
two functions Q v (r) and V^(r) are linearly dependent, that is v v (r) = 
a Q^(r). This means that the function Q^(r), besides satisfying the 
Sommerfeld radiation condition, also satisfies the condition of regu- 
larity at r = 0. It follows that v is an eigenvalue and Q v (r) is an eigen- 
f unction of the following problem: 

Find the values of v for which V(r), the solution of the differ- 
ential equation: 



d 2 (rV) 
dr 2 



[ k 2 (r) - v (v + 1) j (rV) . (4-4) 



is regular at r = and satisfies the radiation condition at infinity, 
that is 



lim I 3- (rV) - ik (rV)| = . (4.5) 
dr o 

It will be seen in later sections that the zeros of the Wronskian, 
v, will be very large, of the order of magnitude of ka. This allows us 
to use asymptotic formulas for the Legendre functions P^(cos (tt - 0)). 
We have 



1/2 



V COS 6) ~ W) sin > COS [(V + 2> 6 " 4 ] (4 " 6) 



for large v, so that 



P (-cos 0) P (cos (tt - 0)) 
v v 



sin vtt sin ttv 



" ( ,(v + 2l l) sin 9 )1/2 ex P { " i ^ + b 6} - 
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Next, using the fact that the dielectric constant of the earth 
is constant, it has been shown [12] that the boundary condition at 
the surface of the earth can be written as 



T_l (rVs ) - r (rVs )l = o 

L -> r=a 



(4.8) 



2 2 

r = i(k/ - k/) 

I O 



1/2 



where k2 and k Q are the wave-numbers of the earth and free space, res- 
pectively. A perfectly conducting earth implies ^ = 00 and hence 
(4.8) reduces to 



(rV s ) 



= 



(4.10) 



r=a 



Equation (4.8) can also be written in a preferred form: 



\ (rV s ) / 



- r . 



(4.11) 



r=a 



We shall use this result presently. 

It is shown in Appendix C, that with the aid of equations (4.10) 
and (4.11), the expression (4.1) can be written for an inhomogeneous 
earth and a perfectly conducting earth, respectively, as 



U(r,9) 



-7^ 

4a v 



(2v+D V b) V r) V cos (1T - 9)) 



sin vtt Q (a) Q (a) 



(4.12) 



dt 



d7 (r Qt> 



r=a 
t=v 
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(2v+l) ( b V b))(rQ v (r)) V COS (ir - 6)) 
sin V7T 



37 <'V] [h 1^)] , 



r=a 
t=v 



Using the approximation (A. 7), we obtain 



(A. 13) 



and 



1/2 



a v 



Q v (a) Q v (a) 



exp ti(v + j) 6] 



3t rQ 



J r=a 
t=v 



(4.14) 



El/2 



U(i 



it sin 6 



Q v (b) Q v (r) 



J r=a L J r= 



a 
t=v 



(4.15) 



Equation (4.14) and (4.15) are the expressions for the Hertz potential 
satisfying the boundary condition (4.11) and (4.10), respectively. 

The formulation in (4.14) (also (4.15)) has the advantage that 
the same formula for U(r,6) represents r > b and r < b, in contrast 
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to eq. (4.1). The expression, however, breaks down for = since 
for small 0, 



P (-cos G) -> log 

whereas the function P^ represents is to be regular for 0=0 and 
r ^ b. It implies that the whole "ray" 0=0 must be considered a 
singularity of the representation in (4.14) and (4.15). This condition, 
however, presents no problem in this investigation, since we will 
always be interested in fields well removed from the 0=0 line. 
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SECTION 5 



THE LANGER APPROXIMATE FORMULAS 



In Section 4 it was shown that the problem of determining the 
Hertz vector U(r,0) and then the electromagnetic field depends upon 
the solution of the eigenvalue problem for the differential equation 
(4.4). In this section we will discuss the "extended" WKB approxi- 
mations due to Langer to obtain solutions of this differential equa- 
tion. The extended form is valid in the neighborhood of the transition 
point as well as at the transition point itself. 

Langer has demonstrated that the solutions of the so called re- 
lated second order ordinary differential equation for a single trans- 

2 

ition point, i.e., a zero of [<J> (z) - w(z)] in eq . (5.1), 

g(z) + [4> 2 (z) - o)(z)] g(z) = (5.1) 



is given by 



(z) = nz) S 1/3 Hf^ (?) j - 1, 2 (5.2) 



1/3 



in which the symbols H^ stand for the Hankel functions of order 



th 

1/3, of the j kind, and 



y(z) = S 1/6 (z) 4> 1/2 (z) (5.3) 



C(z) =/ <|>(z) dz (5.4) 



( H"'(z) _ 

u(z) = "7(zT (5 ' 5) 

2 

where z^ is a zero of [4> (z) - w(z)]. 
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Although the formula for ¥(z) seems to assign a singularity to this 
function at z^, it is found that this singularity is removable, and 
the ¥(zi) is then different from zero. The function to(z) is thus 
bounded in a region about z^. 

2 

If it can be shown that u>(z) is very small compared to <J> (z) 
over the considered region of z, then (5.2) may be considered as a 
solution of eq. (5.1) to practical accuracy. In that case 



g(r) = rQ v (r) , (5.6) 



and 



■ 2 > * , 2 [2. , v(v+l)1 , K 7N 

(f> (r) = k Q k (r) ^ — 2 ' ^ 5 ' 7 ^ 

L k r 



In the region of the zero, if 



U 2 (r)| « 1 » (5.8) 



then the condition 

|u)(z)| « 1 (5.9) 

must be satisfied for (5.6) to be an approximate solution to the 
wave equation. If these conditions are satisfied in the interval of 
interest, for which there is only a single transition point at r = 
then the general solution is given by 



where 



rQ = . 



C =/ ♦ dr (5.11) 
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and <f)(r) is given by (5.7). Thus eq. (5.10) represents the approxi- 
mate solution in the region that includes one transition point. For 
the regions of the second and third transition point, solutions of a 
similar form apply, each valid only in a region including one transi- 
tion point. 

2 

It has been shown that if <f> is linearly varying, oj(z..) = and 

2 

hence eq. (5.9) is satisfied as well as the condition |to(z)| << $ (z) 
for all z. 

Furthermore, Langer has shown that for large Z the asymptotic 
forms of the Hankel functions may be used. Hence 

(1) e H 
H l/3 ™ = (5.12) 

(2) 

H l/3 " C 2^T72 (5 - 13 > 

where and are constants that depend on the normalization of 
these functions. If (5.12) and (5.13) are substituted in (5.10), we 
obtain r 



1 f * dr e 1 f 4> dr 

Jr 1 Jr 1 



rQ = C 7 i + C A -± . (5.14) 



It may be recognized that (5.14) are the "ordinary WKB I? solutions. 

2 

This form of the solution, however, breaks down at a zero of <f) . 

2 

For the discussion of conditions on the permittivity when <f) 
is not varying linearly the reader is referred to Doviak and 
Goldhirsh ' fH] and Bremmer [15], 
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SECTION 6 



SOLUTIONS FOR THREE TRANSITION POINTS 



As we have mentioned in the previous section the extended WKB 
solutions are valid in a region including only one transition point. 
It was pointed out in the Introduction that for the problem of the 
elevated duct, three transition points exist. Hence we must obtain 
three separate solutions for each of the three regions. These solu- 
tions are then matched at points between the transition points using 
the continuity boundary conditions. 

We will represent the solutions (rQ^) to the homogeneous wave 

equation in regions I, II and III by ( r Q v )^» ^ r %^2 aTU * ^ r ^v^3* res " 
pectively. Region I extends from r = a to r = r\^, region II for 

n-, < r < n » and region III from r = ru to infinity. Since (rQ )3 
± - - z z v (2) 

must satisfy the Sommerfeld radiation condition, only solutions 



are allowed. Hence, we have in Region I: 

1/2 1/2 



in region II; 



5 1/2 c 1/2 



in region III: 



1/2 

I 

'3 



W»>3 - F v <T> «ul «3> ■• (6 - 3) 



where 



r 

C x (r) - I <f>(r) dr 
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C 2 (r) = / <fr(r) dr (6.5) 



5 (r) =/ dr (6-6) 

r 3 

f 

(rQ ). and (rQ ) . (j = 1,2,3) (prime means differentiation with res- 
pect to r) must be continuous at r = and r = n . The coet f ici «nts 
, B , and may be expressed only in terms of F^. Further. 1 «- r 

1/2 

,(l) ,^lnls „(l) v 7 . 

*lni = ( "*7 } H l/3 (5 lnl } ffi ' 7 ^ 

1/2 

(2) = /^nl. (2) ( . ( 

A lnl ( * t } 1/3 u lnl ; 1 ' 

1/2 

; 1/2 

1/2 

(2) _ , 2n2 . (2) , . (6 

X 2n2 " ( <f> 22 } "l/3 (C 2n2 ) (b " 1 - ) 

1/2 

^2- <<^> H $ ( ^3n2> < 6 '"> 



where 



r = I 4>dr (6. U) 

lnl 
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r 2 

€ m 
I 



,n 2 



n 2 



? 3n2 =/ * dr • (6.17) 



Using the above in equations (6.1) to (6.3) and applying continuity 
conditions at and we have the set of four simultaneous equations 

A \^}-\ + B A n (2) = C xj 1 } + D X„ (2) (6.18) 
v lnl v lnl v 2nl v 2nl 



v lnl v lnl v 2nl v 2nl 



(6.19) 



C X„ (1 J + D A^ 2 } = F x[ 2 \ (6.20) 
v 2n2 v 2n2 v 3n2 

C X 9 (1 >' + D X 9 (2 >' = F X< 2 9 ' (6.21) 
v 2n2 v 2n2 v 3n2 



which gives for A and B 
° v v 



A = K F (6.22) 
v v i v 

B v *\>2 F v (6.23) 
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where 



+K 



r x (2) 



(2) 

l 3n2 



L, + X 



(2), ,(2), 
lnl 3 n 2 



(2), ,(2) 
lnl 



,(2) (2), 
lnl 3n2 . 



J 

(6.24) 



where 



D - 



.(1) (2), .(1), .(2) (1) (2), (1), (2), 

X lnl X lnl " X lnl "lnl X 2 n 2 X 2 n 2 " A 2 n 2 X 2 n 2j (6 * 25) 



and 



A 2 l h + >i / ^ s^ 1 }' A f 2 s^! 1 ! ^* L . 

lnl 3n2 1 lnl 3n2 2 lnl 3n2 i lnl ^ 4 



where 



D 



(6.26) 



, (1). ,(2), _ , (2), (1), 
v 2 n l 2 n 2 A 2 n l 2 n 2 



(6.27) 



.(1) (2) .(1) .(1) 
A 2 n l 2 n 2 A 2 n l 2n2 



(6.28) 



. (2) (1), _ (1) (2), 
A 2nl 2n2 A 2 n l 2 n 2 



(6.29) 



.(2), (1) _ (1), (2) 
L 4 ~ A 2nl 2n2 A 2 n l 2 n 2 



-..30) 



The functions L , l^, and L 4 de P end onl y on the l a Y rr between 
and It will be seen later that this interval < i" <- n., also 

corresponds to the tropospheric layer. 

Substituting (6.24) and (6.26) into (6.1) gives 
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(rQ ) = F 

v v 



z 1/2 ? 1/2 

\ <^> H V3 + K v 2 ^ «{/] «=l> 



(6.31) 

Substituting this in the boundary condition at the earth, we obtain 



_9 
9r 



♦l 



1/2 



1/2 



H i/3 "l> " k ( ^ 



c 1/2 c 1/2 

^ - k <^> H V3 «1> 



= r 



r=a 



(6.32) 



where 



k = -K /K . 



(6.33) 



In the case of a perfectly conducting earth T = «>, and equation 
(9-32) reduces to 



1/2 1/2 

H $ «i> - k 4)1 <V 



= o 



(6.34) 



r=a 



or 



H l/3 "l> 
4/3 < 5 1> 



= k 



r=a 



(6.35) 



If (6.22) and (6.23) are substituted into (4.14) the coefficient F 
in no way enter the expression U(r,8). 



SECTION 7 

MODIFICATION OF THE RADIAL 
DIFFERENTIAL EQUATION 

The next step is to obtain solutions of the homogenous differential 
equation (3.5). We will write equ. (3.5) in the following form: 



d^u du /.in 

n , 2 n, ri 2 /N n (n+l) 1 n 

-^- + - [k (r) - — ] u n = . (,.1, 

dr r 



Letting h = r-a and substituting in (7.1), we have 



a 2 2 , 2 

d u „ „ a k 

— ^ + [k Z \T (h) * ] u = yl.2) 

dh 2 ° (h+a) 2 n 



2 2 

where we have written a k for n(n+l) 

n 



Now 



a 2 k 2 

^^k 2 (1 , (7.3) 



(h+a) 2 ' n 



so that equation (7.2) becomes 



d 2 u 



dh 



^ + [k 2 k 2 (h) - k 2 + 2k 2 -] u = 
2 o n nan 



(7. A) 
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Since usually, 

k(h) « 1 

and 



k = k 
n o 



the part in brackets in equation (7.4) becomes 



k 2 k 2 (h) + 2k 2 - - k 2 N 2 (r) 
o n a 



2 2 

to within terms of the order of h /a , and 



Letting 



and 



N(r) - Igtf. 

a 



N 2 (h) = 1 + y(h) 



y(h) = 2 x l(f 6 [M(h) - M ] 



2 2 
k - k (1 - X ) , 
no n 9 



we obtain for equation (7.1) 



,2 

d u 9 

— + k Z [y(h) + X ] u « . 
.,2 o J K / n J n 

dh 
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The effect of this transformation has been to transform the 
spherical earth problem into the so called "flat earth" problem. The 
resulting differential equation (7 ,12) turns out to be more easily 
solvable. In the case when the index of refraction varies linearly 
with height, the solutions of equation (7.12) take on particularly 
simple forms. The price one pays is that the solutions will have an 
error of approximately 2% for the frequencies and elevations of interest 
here. This we will tolerate. 

In terms of the function in the general solution (5.10), equ. 
(7.12) becomes 

d 2 Q (h) 

\-+ k Z [y(h) + X 1 Q (h) = . (7.13) 

dh 2 V V 



The general solution of equ. (7.13) is 

1/2 



where 



r h 

? (h) = / 4>(h) dh (7.15) 
J h. 



and 



♦(h) = k Q [y(h) + X v ] . (7.16) 

Next, we will write the Hankel functions of order one-third, H^^(^) 
in terms of the Furry modified Hankel functions h(z) . Some of the 
properties of these functions are given in Appendix E. Thus, the 
transformation is 
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h i <z) = ? 1/3 (h) (5(b)) 



(7.17) 



where 



2/3 



z = (f 5(h)) 



(7.18) 



3/2 



C(h) = j z 



(7.19) 



and equ. (7.14) becomes 



V h)= T# [A v h l (Z)+B v h 2 (Z)] 
<P (h) 



(7.20) 



This form of the solution to equation (7.1) will be used in all sub- 
sequent calculations. The remaining equations of Section 5 and Sec- 
tion 6 are transformed in a straight forward if somewhat tedious way. 
In particular equation (6.35) now takes the form 



h 2 (z(h)) 
h 1 (z(h)) 



- k . 



(7.21) 



h=0 
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SECTION 8 



MODEL FOR THE TROPOSPHERE 



As we have seen previously, a smoothly varying index of refraction 
with an elevated layer will, in general, result in the occurrence of 
three transition points. Instead of an analytic function for the index 
of refraction, we will choose three linear line segments which inter- 
sect at the heights h^ and as in Figure 8.1. The heights h^ and 
h^ correspond to those of n-^ and i n Section 6. These are the 
height at which we "joined" the solutions in the three regions. 

is the height of the layer and (h^ - h^) is its thickness. 
Mathematically, the model can be described as, 



y(h) - 



y Q + a x h 



1 h 1 \ 



y L + a 2 (h - h 1 ) h 1 < h < 



y 2 + a x (h - h 2 ) h ^ h £ 



(8.1) 



We have chosen the same slope (standard) for y(h) in region I and III, 
The slope a^ in region II is simply given by 



3 2 = 



y(h 2 ) - yG^) 
h 2 " h l 



(8.2) 



The model y(h) is illustrated in Figure 8.2. 

For this choice of model the variables t, in equ. (6.4) to (6.6) 
and (6.14) to (6.17) take on the form 



Mq 

Figure 8.1 LAYER MODEL 



M(h) 




and 



:i (h) = k o y* 



1/2 

[y Q + a^h + XJ dh (8.3) 

h 



01 

1/2 



C 2 (h) = k Q I [y x + a 2 (h - h 1 ) + Xj dh (8.4) 



C,(h) = k / [y + a (h - h ) + X ] dh, (8.5) 
J o / l 1 z n 

h 03 



1/2 



? k 3/2 
C ll (h l> = 3 <^> [ *1 + X n ] < 8 " 6 > 

k 3/2 
^21 ( V " 3 <^> ^1 + X n ] < 8 ' 7 > 

, k 3/2 
C 22< h 2> "I <if [ ^2 + X J ( 8 - 8 > 

C 32<V = ! <^> ^2 + X n ]3/2 - ( 8 - 9 > 

In terms of the variable z defined in equ. (7.18), we have for equ. (8.6) 
to (8.9) 

k 2/3 

Z ll ' & ^1 + V (8 - 10) 
k 2/3 

Z 21 = ^ <*1 + V (8 - U) 
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k 2/3 

Z 22 " <4> <?2 + V 
k 2/3 

z 32 = <^> cy 2 + V 



(8.12) 



(8.13) 



Since, for the model of Figure 8.2, is negative, we will choose the 
principal root, that is 



-1*L k 2 /3 



'21 



(8.14) 



and 



-i 2JL k 2/3 

..... 3 1^1 Cy 2 + * n > 



'22 



(8.15) 



In terms of this model, the solution Q (h) , equ. (7.20), takes the 
following form in region I '• 



1/6 



Q v (h) = 

3a l k o 



2/3 

Vl (<^> ( ?o + a l h + v) 



/k 2/3 



(8.16) 



for < h < h 



1 ' 



in region II : 
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2 o 



C v h l <*i + a 2 (h - V + v) + 



(8.17) 



for h 1 £ h £ h 2 ; 



and in region III : 



Q v (h) 



1/6 



2/3 



<— TT> F v h 2 ( ^ (y 2 + a i (h 

3a.k \ 1 

1 o 



h 2 ) + x v )i 



(8.18) 



for h >_ h 2 . 



Expressions (8.16) to (8.18) are substituted in the transcendental 
equation (7.21) and solved for the eigenvalues X . 
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SECTION 9 



EXPRESSION FOR THE ELECTRIC FIELD 

In terms of the function and variables defined in Section 7 and 8, 
the expression for the potential U(r,9) in equation (4.13) takes the 
form 



Uf h -M \^ 2v + 1 V n o 

u(h,e) - —2 Z^ii^e ~rr 



Q (h ) Q (h) P (cos(ir - 6)) 



h=0 



(9.1) 



37 V h > 



h=0 
t=v 



where v is related to X by 



V " k o a(1 - 2 V 



(9.2) 



and 



9-2 9 



9t k a 9A 
o v 



(9.3) 



and h refers to the height of the transmitter, 
o 

The coefficients A and B can be related to k using equ. (6.22) 

v v 

and (6.23) and (6.32). Thus the expression for the height function 
Q (h) , equ. (8.16) becomes in region I 



Q (h) = F K ( 

V V 2 3a, k 



1 o 



l/6 r , 2/3 



h 2 ( ^ ( y Q + a l h + V) " 



(9.4) 



/k 2 ' 3 N 

k h i \£> ( y Q + a i h + K) 



for 



< h 



Similarly, the coefficients C and D can be written in terras of 

v v 

F K . However, in the present investigation we will limit our- 

v V2 

selves to determining the fields above and below the layer. For 
h _> h^, we have 




(9.5) 



for h >_ h 2 



Substituting equation (9.4) into (9.1) and carrying out the necessary 
differentiations, we obtain the following expression for the potential 
below the layer, 



-m v v + i 

, 2 4/3,. vl/3 sin vtt v 

4a a (k p) v 
o 



(h 1 (0)h 2 (h Q ) - h 2 (0)h 1 (h Q ))(h 1 (0)h 2 (h) - h 2 (0)h x (h))/ 



1 - 



a(k ap) 
o 



H>73- k ' h i 2 (0) 



(9.6) 



for 



<_ h < h 



45 



The notation h^CO), h 2 ^ h 0^ and h i^ h ^ refers to the value of the Hankel 
function at the values of h = 0, h = tig and h = h, respectively, in equ. 
(9.4). In this derivation use was made of the Wronskian 



h l (? l } V (? 1 } " V (? 1 } h 2 (? 1 } = " ai (9 ' 7) 



where a = 1.457495. The coefficient p is related to a^ by 



1 pa 



p is equal to (2/3) for the standard (4/3) earth atmosphere. In equ. (9.6) 
k 1 is the derivative of k with respect to the eigenvalue A^. 

Since v is of the order of k Q a, we can make use of the asymptotic 
expressions for the Legendre functions of large order. 

In that case, 



i 1/2 

(v + ■=■) P (cos (ir - 6)) 2k a 

j = ( j q) exp 

sin vir ir sin 6 



if 1 -! (v + j) 6 



(9.9) 



In future discussions it will be of interest to compare the potential 
U to its free space value. The value of U is 

w -ik(b)R 

U =|- (9.10) 

o 4tt R 

If we approximate the distance R by a6 and make use of equ. (9*9) in equ. 
(9.6), we obtain 
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U ^sin 6 
o 



1/2 



3tt 



- 1 



) 



2tt 



(k a) 
o 



1/3, 



1/2 



2 1/3 
a p 



i ^ k a 
2 o 



e x 



(h 1 (0)h 2 (h Q ) - h 2 (0)h 1 (h Q ))(h 1 (0)h 2 (h) - h 2 (0)h 1 (h)) 



1 - rk 1 ,2/3 k ' h l 2 < 0) 
a(k ap) 
o 



(9.11) 



for <_ h <_ h^. The eigenvalues X^ are obtained from 



(9.12) 



where k is given by equation (6.32). 

The electric field is obtained from U by using the expression 

9 



(9.13) 



in equ. (2.2). Since 6 appears only in the asymptotic expression (9.9), 
we have 
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_8 
36 



f 2k Q a ' 
ir sin 6 



1/2 



exp 



3ir 



- i (v + |) 8 



(9.14) 



2k a 
o 



\ir sin 6 



1/2 



-ik a 
o 



exp 



3ir 
1 f- 



i (v + p e 



Using equ. (9.14), we obtain for the electric field , 

<f> 



Mk 2/3 <oy / 2k 2 \ V2 i fL - i 1 e 

E o f o I 4 2 

+ " 4a 2 a 1/3 o 1/3 W sln V 



E 



i ^ k a 6 
2 o 



(9.15) 



(h^O)]^ (h Q ) - h 2 (0)h 1 (h Q )) (h 1 (0)h 2 (h) - h 2 (0)h 1 (h)) 



1 k' hl 2 (0) 



a(k ap) 
o 



for <_ h <_ . 

It should be noted that the ratio U/Uq does not differ appreciably from 
the ratio of electric fields E./E. . As long as the distance R at 
which the field is measured is very much greater than a wavelength it 
is sufficient to use the ratio U/Uq. 
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It will be of interest in later calculations to compare the 
value of the potential U in a layered atmosphere with those values 
when no layer is present. For an homogenous troposphere we then have 



„ .1/2 " 1 £ 8 2tt (k a) 1/3 9 1 *L 

U y e o 4 

u " < • Q ^/ 2 1/3 6 

o (sin 9) p 



\ v i v h 2 (h Q )h 2 (h 



e 2 ° 



[y«»] 



and 



(9.16) 



where the are roots of the equation 



k 2/3 

h 2 t-f) X v ) = (9.17) 



/k 2 ' 3 \ 
h 2 (h )= h 2 % + a l h + \f < 9 ' 18 > 



As a check on these results, it is of interest to compare the 
term in the denominator of equation (9.11), namely 



1 - 1 ,2/3 k ' h l 2( V (9 ' 19) 

a (k ap ) 
o 



with the result of Furry [16] for completely trapped modes. Here we 

have written h., (A ) instead of h, (0) to emphasize the fact that X is 
1 v 1 r v 



49 



complex. This factor is proportional to the normalization constant 
for the eigenf unctions , i.e., the solutions of equation (7.13). For 
this factor Furry gets, in terms of our notation 



-2 



-1/2 

(y(h) + a v ) dh - 



-1/3, 



V) Vlwl 2 , 



01 



(9.20) 



where the eigenvalues A are written as 



X = a + ig 
v v v 



(9.21) 



Carrying out the integration, we get 



-1/3 



-2 



2 2 

,P a j 



KM' 



1 - 



2k 1/3 (ili 

° \ a l a 2 / 

, 2 2.1/3 
a(p a ) 



(9.22) 



(yi + a v ) 1/2 /h l( a v ) 2 



We will concentrate on the factor in large parenthesis. For completely 
trapped modes we obtain 



4h /a 9 - a, \ 



3/2 



(9.23) 



_ ix/S £ " \ "1"2 



k ^ e 
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so that expression (9.19) becomes 



1 - 



2h 1/3 (flU* V 3 
a ^ a x a 2 / >1 v 



1/2 



+ 1 3 



a(ap) 



2/3 



(9.24) 



_ ± "o /!2_^A 3/2 2 

3 \ a 1 a 2 / 1 v 1 v 



But at the eigenvalue for completely trapped modes (see Section 10.2) 



2k o f a 2 a l\ , , ,3/2 



3 \ a l a 2 



(yi + X ) - (m - 1/2)tt m = 1,2,3,... (9.25) 



Then expression (9.24) becomes 

1/3 /a - a \ 1/2 

1+ lt> ' W • 



For the completely trapped modes = a^, to a high approximation. 
Furthermore, Furry introduces the correction factor 6 so that 



.tt/6_ . 6 

e h^(a^) = pure imaginary • (9.27) 



Using this fact in equation (9.26) we obtain, 
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2k 



1/3 



1 - 



fa 2 - 
\ a l a 2 y 



1/2 



a(ap) 



2/3 



IV"v>l 



(9.28) 



which agrees with equation (9.22) for the term in brackets. 
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SECTION 10 



NUMERICAL RESULTS 



10.1 Description of the Computer Calculations 

The computations were done in three major steps. First, the 
Hankel functions were calculated; second, the roots of the transcen- 
dental equation (7.21) were found and finally the Hertz potential 
function, essentially a sum of the roots, was obtained. In the ex- 
pressions for the roots of the transcendental equation and the Hertz 
potential or electric field, Hankel functions of the first and second 
kind of order one-third and complex argument z, i.e., H^^^(z), and 
their derivatives, appear. The variable z takes on the whole range of 
complex values, from very small to very large. To evaluate these func- 
tions, we made use of the Furry [16] infinite series and asymptotic 
series expressions for the modified Hankel functions of order one- 
third, h(i)(2) I- see Appendix D and Section 10.2]. As our condition 
on the accuracy of the calculations, we required agreement with the 
tabulated results. This meant that the modified Hankel functions had 
to be calculated at least in double-precision (fourteen place accuracy) 
on the IBM 370-155*. (This was also the limitation of the present 
software available.) It turned out that the agreement of the calcul- 
ated results with the tables was good to seven (out of eight) signifi- 
cant figures after the decimal point. This was satisfactory. The 
series expression for the Hankel functions were used for |z| < 6 and 

the asymptotic expressions for |z| >. 6. The series in the asymptotic 

, -21/2 _ 

expressions was carried to powers in z of z . The agreement between 

the two series at z = 6 is one-to-one. In the case of the asymptotic 

series one must be very careful so as to use the proper expressions in 



As an historical note, it is interesting to reflect that the tabulated 
results were obtained by the Harvard Computation Laboratory (now Aiken 
Computational Center) in 1944. 
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the various domains of the complex plane. In particular, if the 
magnitude of z is large* and z is very near the negative real axis, 
it is more convenient to use the Miller asymptotic formulas valid for 
f < arg <£L [16] . 

The next major task was the calculation of the roots of the 
transcendental equation (7.21). This turned out to be the most diffi- 
cult and time consuming operation of the whole program. The problem 
was complicated not only by the fact that the roots were in the complex 
plane but also by the fact that roots of large real part and extremely 
small imaginary part (trapped modes) had to be obtained. The expression 
for the transcendental equation (7.21) can be written either as 

F " W ~ k(h l>W h l ( V = ° (10 ' 1) 



or as 



h (A ) 

F =wT7- k( V h 2>V = 



(10.2) 



To obtain the roots corresponding to the trapped modes it is essential 
to use the second form for F, since both 1^ and h^ are very large 
simultaneously. The secant and the Newton-Raphson iteration methods 
[17] were used to solve F = for A^. The experience of previous 
investigators attempting to solve similar problems had been that the 
two methods give, in general, good convergence in the complex plane. 
Each method, however, was cited to give no convergence for particular 
geometries of F in the complex plane. The recommendation in those 
instances has been to go simply to another method. It would be pre- 
sumptuous to make a general conclusion on which of the two methods 
gives more rapid convergence, for this depends on the particular pro- 
blems involved; speaking only for our problem it was found that the 
Newton-Rhapson method gives somewhat better rate of convergence. The 
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roots found, however, by each method are identical for the same error 
bound on X , i.e., | *^ n ^ - ^| < e . The disadvantage of both 

methods is the requirement that the initial guess of a root must be 
fairly close to the actual location of that root to obtain convergence 
at all. The advantage of the secant method is that it involves only 
F whereas the Newton-Rhapson method involves the derivative of F, which 
in our problem is not a trivial calculation. Since we were guided by 
minimizing the computer running time, all roots presented in this re- 
port were found by the Newton-Rhapson method. 

The calculation of the Hertz potential presented no special com- 
putational problems, and hence needs no further discussion. A more 
detailed description of the computational work described above, in- 
cluding computer flow charts, will appear in a companion report. 

10. 2 Calculation of the Characteristic Values for the Propagation 
Modes 

The complex roots or characteristic values of the transcendental 
equation (7.21) can be thought of as propagation modes similar to those 
encountered in waveguide theory. However, the boundary value problem 
of the elevated duct is much more complicated than that since the ele- 
vated layer can not be simply modeled as a perfectly reflecting stratum. 
The modes that actually exist can be classified as "trapped 11 modes, 
"transitional 11 modes and "leaky" modes. This becomes clearer if we 
think of the geometric interpretation in Figure 10.1. The first trapped 
mode is m - 1 and -the last is m = 5 . There are always a finite number 
of these. In terms of the roots X^ of equ. (7.21), 

h,(X ) 

1 V 

This means that the imaginary part of X^ is extremely small and to a 
high approximation X^ = -y(h^). Mode number m = 5, for example, is 
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confined essentially to the heights between h^ and h^ and is char- 
acterized by a very small variation of attenuation with range. When the imaginary 
part of A^, Im(A^), is small and the real part Re(A^) ^ -yO^), we are 
said to have transitional modes, e.g., m = 6 in Figure 10.1. Physically 
these modes correspond to waves that just leak out of the duct. The 
characteristic values A^, that have real parts comparable in magnitude 
to their imaginary parts are said to be leaky. These modes do not 
allow such a simple geometric interpretation as the previous two. 

Now let us turn to the solution of equation (7.21) for A^. We 
started by finding the roots corresponding to the trapped modes. But, 
as we pointed out in Section 10.1, for rapid convergence to a root by 
either the secant or the Newton-Rhapson iteration methods, it is 
essential that the initial guess of a root must be fairly close to the 
root itself. To get a handle on this problem we approximated the 
Hankel functions appearing in (7.21) by their asymptotic representa- 
tions, keeping only the first terms. In this case k(* v ) is given 
approximately by 

4k a 9 - a- 3/2 

i 3 12 
k ^ e J e (10.3) 

where the notation has been explained previously. And equation (7.21) 
becomes essentially 

1 + e " iX = (10.4) 



where X is given by 

4k a - a 3/2 

x = -r [ *i + V • (10 - 5) 

Equation (10.4) is satisfied for 
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X = mTT (m is odd). 



(10.6) 



As a check, we compared this result with the Furry -Gamow phase- 
>ral 
mined by 



integral method [ 4 ] . In this method the real part of X^ is deter- 



h 02 1/2 
k a / [y(h) + X v ] dh = (m - \) tt (10.7) 



, m = 1,2,3. .. . 

01 

where the h^ and h^ satisfy y(h) + X^ = 0. A straight-forward 
integration in two parts from h^ to h^ and h^ to h^ gives 



2k a a 2 " a l 3/2 1 

( I - L ) ly, + A] = (m - h tt (10.8) 



3 a 1 a 2 J l 



m = 1 ,2 , 3 . • • 



It is seen that the results are identical. Equation (10.7) is known 
as the "phase-integral condition 11 . The cases for which this method 
applies are as follows: 

(i) The curve y(h) must have at least one minimum, in our 
case it has two, y(h^) and y(h^)y a mode to which the 
method can be applied must have 

-yO^) < Re(X v ) < -y(h 2 ) (10.9) 

so that the waves are trapped; and it must have 

Im(X v ) « -(y( hl ) + Re(X v )) aQ1Q) 
ImU ) « -(y(h„) + Re(X )) ; 
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(ii) The value of Re(X^) is found from equ. (10.8). 

After the Re(X ) is known, Im(X ) is calculated 

v v 

from an explicit real expression to be derived 
presently. 



For the determination of the imaginary part of X^ we will use 
the modified index of refraction profile y(h) of equ. (8.1), Figure 

8.2. Let us write the characteristic values X as 

v 



X = a + ig . (10.11) 

V V V 



We will now extend the formulas Furry [3] derived for a single minimum 
of y(h) to t\ 
Figure 10.2) 



of y(h) to two minima. The imaginary part 3^ is then given by (see 



o " 2 w 

B=Cke V (10.12) 
v v a 



where 



h ° 2 1/2 
C v " 2 k 2 / (y(h) + o v ) dh (10.13) 



-/ 
/ 

•/ v, 



01 

h 03 1/2 

W = k / |y(h) - a J dh (10.14) 

m a 



v 1 

h 02 



and k = 2tt/X, where X is the wavelength. When equs . (10.18) and 
a 

(10.14) are applied to the profile of Figure 10.2 we obtain 

-2 1/2 a 2 " a l 

C v - 4 (y x + « v ) (-^-i) (10.15) 
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V 2 Y, Y(h) 



CD 
O 

in 

rO 



Fl gure 10 .2 Y ( h ) PROFILE 
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and 



2k 3/2 a 9 - a- 

w m = ~f ( " y 2 " \> (J 7-^> • < 10 - 16 > 

The are determined from equ. (10.6) or (10.8). For later comparison 
we shall call these a^, B y the Furry-Gamow roots. It is seen from 
equ. (10.16) that the condition (10.10) requires 

-2 W 
m 

e « 1 

and (10.17) 

2 

C « 1 . 

v 



This implies that the values of 8 will be in error near a = -y., and 

v v J l 

= -y . It is surprising, however, how good these approximate values 

of a and g are to the exact values determined by the Newton-Rhapson 

v v 

method. As a particular example we have chosen the following constraints 



X = 1 meter 

yi = 211.66 x 10~ 6 

y = 131.66 x 10~ 6 



a. = 2.3518 x 10 7 



a„ = -8.0 x 10 7 



In Table 10.1 are shown the Furry-Gamow roots and their exact values. 
From these results it is seen that the real part of X is in error by 
less than one percent. The imaginary part of X^ is in error by about 
10%; except for m = 2 where the error is about 30%. The 10% error can 
be explained from the fact that condition (10.17) is not always satis- 
fied. There is no ready explanation for the 30% error in mode number 2. 
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It can be concluded, however, that the Furry-Gamow values serve as an 
excellent guide for the initial guess of a root and for this reason 
alone are indispensible . 

There are no such relatively simple formulas that give the roots 
corresponding to the transitional and leaky modes. It was fortuitous 
that the location of these roots could be predicted reasonably accur- 
ately by linear extrapolation from the trapped mode root values. How- 
ever, in only one instance could the location of three successive roots 
be predicted by linear extrapolation. The reason for this failure can 

be seen from the non-linear character of 3 as a function of mode 

v 

number in Figure 10.3. It is perhaps advantages to use some curve 
fitting routine to predict more root values in advance; this was not 
tried here. Such a procedure would certainly cut down on the computer 
running time if it proves successful. 

There was another way in which we went about locating a root 
for a particular set of values of the parameters A, y^, y^, a^ and • 
This approach was motivated by Northover's [18l investigation in 
which he considered reflections from an elevated layer by treating 
such a layer as a discontinuity in the index of refraction. He pre- 
dicted that for the important modes in ground to ground communications, 
namely those waves that are reflected by the elevated layer back to 
the receiver should have corresponding characteristic values whose 
magnitude is less than unity. In terms of our notation this condition 
means that 

v2/3 

< 1 • (10.18) 

They correspond to leaky modes. It turned out that these roots could 
be found rather quickly since there were always just a few of them 
(see e.g., Table 10.2). And the real and imaginary parts were of com- 
parable magnitude. These roots then provided another starting point 
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for finding additional characteristic values. This procedure proved 
to be very convenient. 

10. 3 Discussion of the Characteristic Values and Their Physical 
Interpretation 

All the characteristic values. X , we have found are listed in 

v 

Tables 10.2, 10.3 and 10.4. The accuracy of these values obtained on 
the computer was actually 14 places; they have been rounded off to five 
places after the decimal in the tables. The values were obtained as 
a function of frequency (wavelength X), and layer height, h_^ , as a 
function of lapse rate within the layer, a2 and outside the layer, 
a^ and as a function of the value of y(h) at the start of the layer, 
y^ and at its termination y^ . 

The real part of the characteristic values enters as a phase 
factor in the final expression for the field, equ. (9.15). Each 
characteristic value or mode is seen to have a different phase and 
hence, the addition of many such modes results in appropriate phase 
vector addition. The imaginary part of X^ corresponds to the attenua- 
tion of the mode. This can be written as 

1/2 6 k d 
v o 



attenuation = 




(10.19) 



It is seen from the tables that the trapped modes have a very small 
imaginary part so that the attenuation of these is like (l/yj~d). It 
is noted from a comparison of Tables 10.2 and 10.3 that the thicker 
layer (100 meters) has the smaller imaginary part for the first mode, 
hence the attenuation of this mode is smaller. It will be seen in 
Section 10.4 that this will also correspond to more energy being 
trapped in the duct. In both cases the first five modes are trapped. 
In Table 10.4, the layer is even thicker (150 m) but we have lowered 
the frequency to 150 MHz and as a result the first mode has an imaginary 
part several orders of magnitude greater than for 300 MHz. In addition, 
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Mode 










Number 


Characteristic 


Values 






m 


Real Part (x 10" 5 ) 


Imaginary Part 


1 


-19.25747 






io -13 


2 


-17.33864 






10" 10 


3 


-15 .78268 


• / § j -j \j 


"V" 




4 


-14.43457 


.17336 




10 - 6 


5 


-13.20194 


. 84462 


1*; 


10 _ 6 


6 


-12.00625 


.16483 


x 


IO" 5 


7 


-10.85324 


.22246 


x 


10-5 


8 


- 9.75080 


. 26214 


x 


10"5 


9 


- 8.69502 


.29129 


X 


io-5 


10 


- 7.68003 


.31404 


X 


IO -5 


11 


- 6.70058 


.33267 


X 


10-5 


12 


- 5.75236 


.34849 


X 


10-5 


13 


- 4.83184 


.36233 


X 


10-5 


14 


- 3.93611 


.37476 


X 


10 * 


15 


- 3.06264 


.38639 


X 


10-5 


16 


- 2.21008 


.42040 




10-5 


17 


- 1.36162 


.42117 


X 


10 Z 5 


18 


- .49410 


.46243 


X 


10 3 


19 


.43138 


.52265 


X 


10 = 5 


20 


1.43287 


.59324 


X 


10 


21 


2.51086 


.66949 


X 


10-5 


22 


3.66053 


.75062 


X 


10-5 


23 


4.87641 


.83685 


X 


10-5 


24 


6.15302 


.92792 


X 


10-5 


25 


7.48485 


1.02197 


X 


IO" 5 


26 


8.86721 


1.11422 


X 


10-5 


27 


10.29907 


1.19673 


X 


IO" 5 


28 


11.78600 


1.26204 


X 


IO" 5 


29 


13.33784 


1.30910 


X 


10-5 


30 


14.96139 


1.34344 


X 


10-5 



TABLE 10.2 CHARACTERISTIC VALUES 
(A = lm, AM = 40 in 30m, h. = 900m) 
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Mode 
Number 


1 

Characteristic Values 

■ 

i 


m 


Real Part (x 10 5 ) 


Imaginary 


Part 


± 


-19.58371 


.12851 


X 


10-15 


Z 


-17.74367 


.77288 


X 


10-12 


3 


-16.28801 


O O O 1 "7 


X 


10-9 


/, 


-15.06395 


o o c on 

. 22680 


X 


10-7 


c 
D 


-14.00183 


.32445 


X 


10" 6 





-12.98671 


.12363 


X 


10-5 


"7 


-11.94416 




X 


10-5 


Q 

o 


-10.90414 


.30204 


X 
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.49691 


X 


10-5 


15 


- 4.50455 


.48825 


X 


10-5 


16 


- 3.66337 


.48082 


X 


10-5 


17 
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10° 
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26 
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X 
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27 
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X 


10-5 


28 


8.40164 


1.08675 


X 


10-5 



TABLE 10.3 CHARACTERISTIC VALUES 
(A = lm, AM = 40 in 100m, h, = 900m) 
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only the first three modes are trapped resulting in lower signal 
strength within the duct. 

10.4 Comparison of Calculated and Experimental Results 

The calculations of field intensity relative to free-space were 
performed using expressions (9.11) and (9.12). The results are plotted 
in Figures 10.4 to 10.8. All curves are for a frequency of 300 MHz, 
and a layer height of 900 m. Two values of layer thickness, 30 m and 
100 m, but of the same decrement in modified index of refraction, 
AM = 40, have been used. Such layer values are representative of those 
measured in the trade winds ducts near San Diego, Calif. [19] and the 
South Atlantic [201 . The field intensity is plotted as a function of 
height and distance for various transmitter locations. A comparison 
of the plots shows that the thicker layer (100 meters) has also the 
greater field intensity values; a difference of about 10 dB on the 
average exists for the two layer thickness values. This appears logi- 
cal since, for the same decrement (AM = 40) the 100 meter layer results 
in a wider duct than the 30 meter layer and, hence, is capable of 
trapping more energy. When the transmitter is moved out of the duct, 
Figure 10.6, the field intensity decreases by approximately 5 dB on 
the average. In addition, we notice in this case interference between 
the waves that move above and are coupled into the layer and those that 
go through the duct strike the earth and then escape into the space 
above the layer. In Figure 10.6, constructive and destructive inter- 
ference are observed at approximately 450 km and 500 km, respectively. 
The field intensity is given as a function of height in Figures 10.7 
and 10.8. Here again, the wider duct results in a higher field strength, 
about 10 dB. We also notice in Figure 10.7 that we obtain interference 
type (multipath) lobing; the intensity of which decreases as one gets 
closer to the bottom of the layer. Furthermore, the field strength in- 
creases by approximately 5 dB from near the bottom of the duct to the 
top. This latter behavior does not appear to be in evidence for the 
wider duct in Figure 10.8. 
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Figure 10 7 FIELD INTENSITY AS A FUNCTION OF HEIGHT 

(Am = 40 IN iOm ) 
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We will now compare these results with some experimental evidence. 
For comparison, we use the test results of Project Neptune [21]. This 
project involved air-to-air and air-to-ground propagation testing in 
the trade wind ducts at 220 MHz. Its purpose was to explore the communi- 
cation potentialities of the trade winds inversion layers. For compari- 
son we will choose the thinner layer (30 m) . It is to be noted that in 
these test results, the measured values of refractive index values are 
open to question |21 |. The agreement between our results and the tests 
is good on a relative basis. When the results are compared on an absolute 
basis our results are approximately 10 - 15 dB higher. To partially ex- 
plain this difference, we note that our results are for a somewhat higher 
frequency, 300 MHz versus 220 MHz. In general, the higher the frequency, 
the greater the extent of trapping, and hence, the higher signal level. 
Another difference is due to the fact that in our model the duct width 
and AM scenarios constant whereas in the tests the width and decrement 
actually changed in the distance. This would also increase the leakage 
from the duct and result in a lower signal level. In the tests it was 
observed that signal level increased from 5 to 20 dB as the receiver 
moved closer to the bottom of the layer. The field being highebt at 
the base of the inversion. These results agree essentially with ours. 
We obtain an increase in signal level of 5 to 10 dB. Furthermore, it 
was reported that the signal level remained "relatively high" above the 
duct. This also agrees with our calculated results. It being that the 
signal is only approximately 5 dB less above the duct. Finally , the 
tests established good signal correlation for 60 feet ( % 18 m) verti- 
cally. If this result is compared with Figure 10.7, it is seen that 
it agrees quite well with the calculated numbers; our results showing 
about 80 feet at the 3 dB points. 
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SECTION 11 



SUMMARY OF CURRENT RESULTS AND RECOMMENDATIONS FOR FUTURE WORK 

We will here summarize the results we have obtained and explain 
their significance. Also, we will make several recommendations for 
future work. In addition, we will point out the most significant 
elements to be incorporated in a propagation experiment. 

We have shown that the presence of elevated ducts can produce 
significantly higher signal strengths; on the order of 20 dB above 
free space. We have further demonstrated that the tri-linear model 
for the index of refraction is plausible and gives good agreement with 
experimental results on a relative basis but tends to make rather 
high estimates of signal strength on an absolute scale. We have 
shown the existence of interference type lobing at the duct height 
with 5 dB to 10 dB fades at 300 MHz and good signal correlation for 
vertical distances of roughly 80 feet. From an operational point of 
view, it is best to fly the transmitter and receiver at the base of 
the elevated layer. At this location signal strength is highest 
and interference type lobing is least. In order to avoid the 10 dB 
fades it is essential that space diversity reception be employed. 

In the present report we have limited our analysis to the detailed 
analysis of a single frequency, 300 MHz. Future work should be ex- 
tended to a comprehensive comparison of several frequencies. Spatial 
correlation distances of the signals should be obtained for such fre- 
quencies in order to obtain estimates of vertical and horizontal 
antenna diversity distances. Furthermore, a spectral analysis should 
be performed to estimate the signalling bandwidth that can be supported 
in the presence of elevated layers. And finally, the calculated re- 
sults that have been obtained for heights within the duct and slightly 
above it should be extended to include the whole range from the ground 
to great heights. 
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From the above discussion it becomes quite clear what the most 
essential elements of any air-to-air propagation experiment should be. 
The most detrimental effect of this channel is interference type 
lobing. Hence, the experiment must, first, establish the signalling 
bandwidth it can support. This could be performed by probing the channel 
with a multi-tone signal consisting of a group of equally spaced fre- 
quencies. Secondly, the experiment must test some antenna diversity 
system in order to reduce the amplitude modulation caused by fading. 
A general experiment incorporating these elements is described in 
Reference 1. 
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APPENDIX A 



JUSTIFICATION FOR THE HERTZ POTENTIAL 



We shall show here that the one dimensional Hertz vector II = 
(0,0, U) is a satisfactory approximation to the solution of Maxwell's 
equations . 

Let us consider II which satisfies the equation 
— — 2 — 

v • v • n - vv • n - yew n = o (a-i) 

and construct the field vectors corresponding to the magnetic type 
(magnetic dipole), then 



D = - iyetoV • II (A-2) 
H = V ' V " if (A-3) 
The vectors D and H satisfy the equations 

V x E = - iwy H + ^ x E (A-4) 

V x I = itoe E (A-5) 

V • D = (A-6) 

V • B = Vy • H . (A-7) 

It is seen that equations (A-4) to (A-7) do not satisfy Maxwell's 
equations exactly. For the atmosphere, however, the gradient of y 
is extremely small so that terms involving it can be neglected. The 
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ratio of the magnitude of the second term to the first term on the 
right hand side of equ. (A-4) can be written as 



x jm . ffi M . (A _ 8) 

2TT y y 



Except quite close to the source, the field will be almost a plane 
wave for which |e|/|h| is approximately yj y/e . Furthermore, |Vy|/y 
is of the order of 10"^, and therefore the ratio of the second term 
to the first term is of the order X • 10"-* , which is always very 
small compared to unity. From equs. (A-5) and (A-7) it follows that 



V ' B „ X_ | Vyl \/ y £ |H| 
i „ ~i 2tt y e i— i 



(A-9) 



and this again is of the order of X • 10 "\ It follows that equs. (A-4) 
to (A-7) differ from MAxwell's equations by a negligible amount, and 
we may regard equations (A-2) and (A-3) as a valid way to obtain the 
electromagnetic field. 

Because the Hertz vector has only its radial component non-zero, 
it follows that [19] 



V(ye) • V • n = and V(ye) • n * . (A-10) 
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APPENDIX B 

CONTOUR INTEGRAL REPRESENTATION FOR U(r,0) 

Following a method due to Watson, we shall show that the series 
expression for U(r,6), (3.9) can be expressed as a contour integral. 
Consider the series 

(2n + 1) f (n) P (cos 0) (B.l) 
n 

where f(n) is an analytic function of n, We will show that S can be 
represented by the following contour integral: 




(B.2) 



where C = C^ + C^ is a contour that starts at 00 - i6 in the complex 
s-plane, goes below the real axis at s = and then above the real axis 
to 00 + i6. Since f (s - and P g _ ( cos " ©)) are analytic functions 

of s, the only singularities of the integral are poles at those values of 
s inside C for which cos tts = 0, that is s = 1/2, 3/2, 5/2, ... 

Equation (B.2) can be written as the sura over C. and C ? . Hence, 




(B.3) 



where H(s) is the integrand in (B.2). In the second integral of (B.3) 
replace s by -s . This transformation reflects the contour C^ about the 
origin into the contour C«, as well as alters the integrand. 
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Hence , equation (B . 3) becomes 



f f 

- -I H(-s)ds + I 



F = -/ H(-s)ds + J H(s)ds (B.A) 
] 3 ' Z \ 



Adding and subtracting / H(s)ds to (B.4) results in 

: 2 



F = / H(s)ds -/ [H(s) + H(-s)] ds • (B.5) 

A + C 3 J C 3 



Using the fact that 



P s - i/ 2 (cos (7T ~ 0)) = P - s - i/2 (coS (7T ~ 0)) (B ' 6) 



and substituting for H(s) we obtain 



■/ * 



i cos (its) f (S - 2> P s - l/2 (cos (lT - 6 » ds 



+ / - — - [f (-s - h - f .(s - b] ds 

I i cos its [ 2 2 J 



(B.7) 



The contour + may be closed by a semicircle of infinite radius in 
the lower half plane as shown in Figure (B.2). Thus (B.7) becomes 
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F =/ 



13 4 



1 COS 7TS 



f (s - \ ) P g _ 1/2 (cos (tt - 6)) 



ds + I, 



(B.8) 



where 1^ is the second integral in (B.7). 

Now substituting the expression for f obtained in the previous 
section, 



f(s - ^) = 



-M 



V 4ttW 



s - 1/2 



% - l/2 (b) V s - l/2 (r) r ± b 



[ V s - l/2 (b > % - l/2 (r) r > b 



(B.9) 



where Q g _ 1/2^ and V s - 1/2^ re P resent analytic functions along 
the positive real axis. Substituting (5.9) into the expression for 
U(r,6) we obtain 



U(r,6)= 



TTi I 



sds s - 1/2 



COS 7TS W 



C l +C 3 +C 4 



s - 1/2 



Q s -l/2< b > V s -l/2 (r) r ^ b 



V s - l/2 (b) % - l/2 (r) r ^ b 



(B.10) 



And if we replace s by v + y and apply the Caudy residue theorem, then 
(B.10) becomes 



m (v + V p , (cos (77 ~ e)) 



<-r> 



s=v 



Q v (b) V v (r) 



V v (b) Q v (r) 



+ I, 



(B.ll) 
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where v represents the zeros of W g in the lower half of the complex 
s-plane • 

It can be shown that if the earth is assumed to be a perfect con- 
ductor, then the integral 1^ vanishes identically. If the earth is not 
a perfect conductor, 1^ will not be zero. However, it can be shown 
that 1^ is small compared to the first term of the series in (B.ll) 
as long as the dielectric constant of the earth is large compared to 
the dielectric constant of the atmosphere. So that finally, 




Q v (b) V v (r) 



r < b 



< 



Q (r) V (b) 



r > b . 



s=v 



(B.12) 
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APPENDIX C 

DERIVATION OF EQUATIONS (4.12) and (4.13) 
From equation (4.4), consider the following two equations: 



dr 



2 < r V + 



2 t(t + 1) 



(rV t ) = 



(CI) 



and 



(rQ ) + 



dr 



2 



k 2 _ v(v + 1) 



(rQ v ) = 



(C.2) 



Here t is any value and is not an eigenvalue. Multiply (C.l) by (rQ y ) 
and (C.2) by (rV^) , subtract and integrate with respect to r from r = a 
to r = 00 . We obtain 



(v - t) (t+v+1) 



>/ Vv 



dr = rV -r^ (rQ ) - (rQ ) (rV ) 
t dr x v x v dr t 



a 

(C.3) 



Substitution of (4.8) for the lower limit in (c.3) results in its 
vanishing. Hence, considering the right side of (C.3) at its upper 
limit as t v results in 



/ 



V Q dr = lim 
t->v 



(rV t )(rQ v )' - (rQ v )(rV t )' 



v(v + 1) - t(t + 1) 



(C4) 



The numerator vanishes since it reduces to the Wronskian W as t -> v. 

v 

This also implies that 



V = a Q 
v v v 



(C.5) 
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as was pointed out previously. Then applying L 1 Hospital's rule to 
(C.4) results in 



/ 



dW 
t 

V v Q v dr = lim -(2vTl) • < C ' 6 > 



Since the Wronskain of any two independent functions is a constant (C.6) 
reduces to 

dW 



/ 



v v Q v dr = - -^rrh 



a 

Making use of (C.5), we have 



dW V (r) 

(2v + 1) / (Q ) dr (C.7) 



dt t=v V r) / V 



And substituting (C.7) in the expression for U(r,9), we have 



-M V V b) V r) P v (GOS (7T " 9)) 

u(r.e) = — L js " TuTw for r - a ' 

v / (Q v (r)) Z dr 



(C8) 

The integral in the denominator is numerically inconvenient. To 
simplify (C.8) consider the expression (C.3) again 




V dr = 
t v 



v(v + 1) - t(t + 1) 



(C.9) 



and using (C.5) as 
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V t -« t Q t • 



then the normalized integral is 



/ 



2 Q t (r) 

<V dr = lim roo 

t-n> t ' 



(rQ ) (rVj 
dr t 



v(v + 1) - t(t + 1) 



(CIO) 



Differentiating numerator and denominator with respect to t, and 
remembering that when t = v, an eigenvalue, the function V^(r) is 
identical with Q t (r). 

Hence , 



(Q ) dr = lira 



t-vv 





4- (rQ ) - (rQ ) 
dr x v v 


h <r Vl 


3 
3t 


v(v + 1) - t(t + 1) 





(C-ll) 



The denominator of (C.ll) is simply -(2t + 1), The numerator is 



lim 
t-n> 



(C.12) 



evaluated at r = a and r = 00 . Now since the right hand side of (4.11) 
is independent of t, we have 



3 3f (r V 



3t (rV t ) 



= 



(C.13) 



r=a 
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or 



If (r V h (r V 



r=a 



= 



(C.14) 



r=a 



or since when t = v the function V (r) is identical with Q (r), then 

v x v 



(r V If h (rV v ) 



- It (r V h (r V - 



(C.15) 



r=a 



Using this in (C.12) we obtain for the right hand side of (C.ll) 



/ 



(Q..) dr - 



-(2v + 1) 



(C.16) 



The numerator vanishes as t-**> . If Q v (a) 4- 0, that is, the earth is 
not a perfect conductor, then we can write 



<T(r) dr = 



-1 



(2v + 1) 



< r V 2 If 



(C.17) 



r=a 
t=v 



If the earth is a perfect conductor then from (4.10) 



(rV v ) 



= 



(C.18) 



r=a 



and remembering that as t-*-v, = we have from (C.12) that 
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for the perfectly conducting earth. 

Using (C.17) and (C.19) in the 
dition (4.11) and (4.10), we obtain 



expression for the boundary con- 
expressions (4.12) and (4.13). 
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APPENDIX D 



THE MODIFIED HANKEL FUNCTIONS AND THEIR PROPERTIES 



The differential equation 
H 2 

^-j + z u = (D-l) 
dz 



is known as Stokes' equation and general solutions can be written in 
terms of Bessel functions of order one-third [22]. The only singularity 
of this equation is an irregular singularity at infinity. Equation (D-l) 
has solutions h^( z ) and \i^{z) which can be written in terms of the two 
Bessel functions of the third kind of order one-third as follows: 

hl (.) . (f S'h 111 „<» (f ,W, (1,-2) 



, , . ( 2 3/2, 1/3 „(2) ,1 3/2, , n _ 

h 2 (z) = (jz ) H£^ (-^ z ) (D-3) 



^1/3 anc * ^1/3 are tr ^-P^ e va l ue d functions of x, with a branch 

point at the origin. However, the product of the two multiple-valued 

2 3/21/3 2 3/2 

functions (— z ) and C "3 z ) is the single valued function 

h_.(z). The Wronskian of h.^ and h^, WO^, h^) has an identically van- 
ishing derivative, so that 

4i 3 1/3 
W(h , h 2 ) - - ~ (f) 



= -1. 457 495 441 040 461 i (D-4) 



a constant. 
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The solution h^(z) and ^(z) are a ^- so related to the Airy integrals 
Ai(z) and Bi(z), defined for real values of z by 

t 3 

Ai(-z) = - / cos («- - zt) dt (D-5) 




Then 



where 



J i (_z) = \l [ e_t3/3 " Zt + sin (3- " zt ) d t] 



(D-6) 



h^z) = k |^Ai(-z) - i Bi(-z)J (D-7) 

h 2 (z) = k*[Ai(-z) + i Bi(-z)J (D-8) 

Ai(z) - h (- z ) + W-z) (D-9) 

Bi(z) = h^-z) - ^ h 2 (-z) (D-10) 

2k 



k = (|) 2/3 (1 - i^) 



or 

k = 1.310 370 697 104 448 - 0.756 542 874 711 451 i 



and 

4- = 0.286 178 560 638 333 + 0.165 225 269 020 841 i 



and the star denotes complex conjugate. 
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Asymptotic Expansions 

_ , . _ , / \ i- 2tt Att 

The asymptotic expansion of h^(z) for - -j— < arg z < — , 



h^(z)~ a z e 



2/3i 



3/2 _ 5_n_i r 
: 12 



1 + 2 (-i) m C z _3m/2 



m=l 



(D-ll) 



c = (9-4) (81-4) ... (9 [2m-l]' -4) 

m ~4m »m . 

2 3 m! 



This expression cannot be used along the ray, arg z 
is the branch-cut for the multiple valued functions, 
coefficient is 



- y~. This ray 
The numerical 



a = 2 1/3 3 1/6 tt 1/2 = 0.853 667 218 838 951. 



An asymptotic expansion valid on the ray, arg z = - rr- is 



2/3i z 3 ^ 2 - — 
h 1 (z)~a z" 1/A e ' 12 |1 (-i)"' C z 

1 m=l m J 



[l + X) (-J 

L m=l 



/2~| 



(D-12) 



-1/4 

+ a z e 



-2/3i z 



3/2 _ IIttI 
12 



L m=l 



(i) m C z" 3m/2 
m 



This expansion holds for - — < arg z < 0; the branch cut for the 
fractional powers of z can be drawn anywhere within the sector 
< arg z < — . 

The existence of two expressions of different forms which represent 
asymptotically the same integral function, h^z), is an example of Stokes 1 
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phenomenon. 

2tt Att 

The region — < arg z < — for (D-ll) coincides with the region 
4tt 2tt 

- — < arg z < - — for (D-12). Thus in this region, we use arg z = a 

in (D-ll) and arg z = a - 2tt in (D-12) with ~ < a < j^. 

The asymptotic expansions for ^^(z) are 



h 2 (z) ~ a z e 



3/2 , 57ri 
-2/3i z + 



12 [i + £ « 

L m=l 



i) m C z" 3m/2 



m 



4tt 2tt 
valid for - — < arg z < — , and 



(D-13) 



o /o • 3/2 , 5Tri 

h , , -1/4 " 2/31 z + TT 



L m=l 



C z 
m 



-3m/ 2 



(D-14) 



9/ , 3/2 llwl 
2/31 z + — 

+ a z e 



1 + £ ("i) m C_ z" 3m/2 



m=l 



m 



valid for < arg z < — . 

Special Properties of h^(z) and hjCz) 

The zeros of h- (z) and those of h- 1 (z) lie on the ray of slope 

2tt 

- — , while those of \i^{z) and of h^^z) lie on the ray of slope 
2tt 

angle — . The zeros of all these functions are simple. 



If 



T = | z 3/2 and z = (f T) 2/3 (D-15) 
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with 



Then 



and 



arg z = — arg T (D-16) 



(T) = T 1/3 h x (z) (D-17) 
(T) = T" 1/3 h 2 (z) (D-18) 



1/3 -^i 

H 2/3 (T) = ( f } 6 3 T V 00 (D_19) 

1/3 + — 
C2"i 2 3 -2/3 

H 2/3 (T) = ( V e T Z/J h 2 '(z) (D-20) 



Figures D-l and D-2 give the regions of validity of the asymptotic 
expansions of h^(z) and h^(z) . 
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REGION OF VALIDITY 
OF EQU. (D-ll) 



REGION OF VALIDITY 
OF EQU. (D-12) 



Figure D-l . REGIONS OF VALIDITY OF THE ASYMPTOTIC EXPANSION FOR h 2 (z) 





REGION OF VALIDITY 
OF EQU. (DM3) 



REGION OF VALIDITY 
OF EQU. (D-14) 



Figure D-2. REGIONS OF VALIDITY OF THE ASYMPTOTIC EXPANSIONS OF h 2 (z) 
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